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Abstract 

We present the plane-symmetric solitonlike solutions of magnetostatic equilibria by solving the 
nonlinear Grad-Shafranov (GS) equation numerically. The solutions have solitonlike and periodic 
structures in the x and y directions, respectively, and z is the direction of plane symmetry. Although 
such solutions are unstable against the numerical iteration, we give the procedure to realize the 
sufficient convergence. Our result provides the definite answer for the existence of the solitonlike 
solutions that was questioned in recent years. The method developed in this paper will make it 
possible to study the axisymmetric solitonlike solutions of the nonlinear GS equation, which could 
model astrophysical jets with knotty structures. 
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I. INTRODUCTION 



The magnetostatic equilibria are of fundamental interest, since they well approximate 
slowly varying magnetically confined plasma configurations. In systems with the helical 
symmetry (i.e. unification of plane symmetry and axisymmetry) , the magnetostatic equilib- 
ria are described by the so-called Grad-Shafranov (GS) equation (see [l] for a review). The 
GS equation is an elliptic equation for the flux function \l/ with a source term depending on 
two functions of \I/ that can be chosen freely. The solutions of the GS equation are often 
used in theoretical studies in the contexts of the astrophysics or the tokamak physics. 

Recently, the existence of an interesting solution of the GS equation was suggested by 
Lapenta [2j. In that paper, the GS equation with a cubic source term (say, the cubic GS 
equation) in the plane-symmetric case was discussed. An analogy between this equation 
and the cubic Schrodinger equation was pointed out, and the real part of the solution of the 
cubic Schrodinger equation was presented as an analytic "solution" to the GS equation. This 
"solution" is periodic in the y direction, and has a solitonlike structure in the x direction. 
Here, z is the direction of the plane-symmetry Unfortunately, an erroneous assumption in 
this analysis was pointed out [3[] and thus the "solution" of [2[] cannot be accepted. How- 
ever, Lapenta performed a numerical simulation adopting his "solution" as the initial state 
and observed that the system is relaxed to a quasi-equilibrium state which maintains the 
solitonlike structure jj] (until the instability becomes relevant). Then, he claimed that the 
solitonlike solution has practical applicability and relevance. The discussion that supports 
the existence of the solitonlike solution was given in [5|. The growth of the instability of 
the solitonlike structure of this system was simulated in order to discuss the collimation and 
_o fastropllyS1Calj e t J. 

In this paper, we present the plane-symmetric solitonlike solutions of the cubic GS equa- 
tion by performing highly accurate numerical calculations. There are several reasons for 
doing so. First, the question on the existence of the solitonlike solutions of the GS equation 
should be answered. Although the simulations in 4|, |6| strongly indicate the existence of the 
solitonlike solutions and explain the gross features of such solutions, strictly speaking, the 
obtained quasi-equilibrium states are not the solutions of the GS equation since they weakly 
depend on time. Besides, even if we ignore the time dependence, the GS equations which 
the obtained quasi-equilibrium states satisfy are not necessarily the cubic GS equation, i.e., 
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the form of the source term should be different in general. Since the possible existence of 
the solitonlike solutions attracts attentions, it is worth studying these solutions by directly 
solving the cubic GS equation. 

Second, the numerical technique for solving the nonlinear GS equation should be devel- 
oped. As we will see later, there is a difficulty in solving the cubic GS equation such that 
the nontrivial solution of this equation is unstable against the numerical iteration, i.e. the 
standard technique to solve elliptic partial differential equations. The development of such 
a technique is important not only in the plane-symmetric case but also in the axisymmet- 
ric case, because the solitonlike solution of the axisymmetric GS equation is expected to 
have interesting astrophysical applications. In fact, the observations of active galactic nuclei 
(AGN) suggest that the astrophysical jets have the magnetic multiple islands [3|. Several 
models for such knotty jets have been proposed, and one possible direction is to model the 

n n n n q 

knotty jets as the magnetostatic equilibria in the comoving frames [2|, m, la, 18|, 19[. As the 
first step to study the axisymmetric nonlinear GS equation, it is useful to begin with the 
simpler plane-symmetric case where the existence of solitonlike solutions is highly likely. In 
this paper, we give a procedure to realize the sufficient convergence and successfully ob- 



tain the numerically unstable solutions. Here, it has to be mentioned t 



rat some isolated 



axisymmetric toroidal Alfven solitons were numerically obtained in [10|, lU[ by combining 
Fourier transformation and the method of the Green's function. Although their method can 
be applicable also for the present cases, our method is somewhat simpler and easier. 

This paper is organized as follows. In the next section, we briefly review the GS equation 
and introduce the cubic GS equation. The asymptotic behavior of the solitonlike solution 
is also studied. In Sec. Ill, we explain the numerical method and estimate the numerical 
errors. In Sec. IV, the numerical results are presented. Some properties of the obtained 
solution are also examined. Sec. V is devoted to summary and discussion. In this paper, 
we adopt the unit where the vacuum permeability /i = 1 because it can be restored by 
dimensional considerations if necessary. 

II. THE GRAD-SHAFRANOV EQUATION 

In this section, we review the Grad-Shafranov (GS) equation for the plane-symmetric 
case and introduce the assumption on the arbitrary functions that leads to the cubic GS 
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equation. Then we explain our requirements on the behavior of the solution. 

We consider static configurations where the magnetic fields B are embedded in an ideal 
plasma with velocity v = 0. The basic equations are Ampere's law and the force balance 
equation: 

J = VxB, (1) 

J x B = Vp, (2) 

with the Gauss law V • B = 0. Here, J is the electric current and p is the pressure of 
the plasma. Note that in the static configurations, the condition of magnetic confinement 
E = —v x B just indicates the absence of the electric fields. 

We introduce the Cartesian coordinates (x, y, z) and assume the plane-symmetry in the 
z direction. Then, the magnetic field can be given by 

B = e z xV^(x,y) + B z (x,y)e z , (3) 

where *&(x, y) is the so-called flux function and e 2 is the unit vector in the z direction. This 
formula satisfies the Gauss law automatically. By Ampere's law (CD), the electric current is 
calculated as 

J = (V 2 ^)e 2 + (VB Z ) x e z . (4) 
Substituting this formula into Eq. (j2J), we find 

[(V* x VB Z ) ■ e z ]e z = Vp+ (V 2 ^) W + B Z VB Z . (5) 

Since the left hand side is the vector in the z direction while the right hand side is the vector 
in the (x, y)-plane, both sides have to be zero: 

(Vtf x VB Z ) ■ e z = 0; (6) 

Vp+ (V 2 ^)Vty + B Z VB Z = 0. (7) 

Eq. (Ej) indicates V\& || VB Z , and then Eq. (JTj) indicates || Vp. Therefore, the contours 
of B z , p, and ^ should coincide, and at least locally B z and p are given by 

B z = /(*), p = g(V), (8) 

where / and g can be chosen arbitrarily as long as they are regular and g > 0. Then, Eq. ((7|) 
is rewritten as 

V 2 * = -g> - ff. (9) 
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This is the GS equation for the plane-symmetric case. 

In this paper, we consider the situations where the GS equation (J5]) is reduced to the 
following form: 

V 2 ^ = (ajj + (3fo 2 ) . (10) 

Here, a and /3 are assumed to be positive without loss of generality. Since the source term 
has the cubic term, we call this equation the cubic GS equation. Eq. (fTPjl can be derived if 
we choose the functions / and g satisfying the relation 

f + 2g = al^+ l -(3l^ + C. (11) 

Here, C is a non-negative constant and it is zero if all physical quantities B and p decay 
at the distant region. There are infinitely possible choices for / and g, since if / = fo and 
g = go satisfy Eq. ( ITTi) . f = fo + X an d g = go — foX ~ X 2 /2 also satisfy this relation, where 
X = x(^) is an arbitrary function. Therefore, one solution ^ of the GS equation ( fTOl) can 
describe many different configurations. 

It is possible to eliminate ao and /3q from Eq. ffTUj) by introducing the new coordinates 

x := a x, y := a y (12) 

and the rescaled function 

u := & ¥. (13) 
By these transformations, the cubic GS equation becomes 

U )S X + Uyy = ~U(1 + U 2 ). (14) 

We require u(x, y) to be periodic in the y direction, to have the mirror symmetry about 
the y axis, and to become zero at x — > oo. Namely, we look for the solutions which behave 
solitonlike in the x direction. 

Let us study the asymptotic behavior of u(x,y) at x — > oo. Since we require that u(x,y) 
decay in this limit, Eq. (IT4T) is approximated as u^x + u^y = —u. A solution to this equation 
satisfying the above requirements is 

u = Aexp {^—x^/q^ 2 — lj sin (q^y) . (15) 

Here, A and q are constants and the value of q is limited as < q < 1. The period 
in the y direction is yp = 2nq. This asymptotic behavior ( |T5l) satisfies u(x, 0) = and 



u } y(x,yp/A) = 0, and we assume that these properties are held for all values of x. These 
two relations together with the condition for the mirror symmetry u t x(0,y) = will become 
the boundary conditions in the numerical calculation. Once the solution satisfying these 
boundary conditions is generated, the solution in the whole region of (x, y) is obtained by 
the relation u(x, y) = —u(x, —y) = u(x, yp/A — y) — u(x, y + yp) = u(—x, y). 

Note that Eq. ( fl5l) is the exact solution of the GS equation f TTDT) in the case ao = 1 and 
(3q = 0. Therefore, without the cubic term in Eq. (Tbll) . the solution diverges at x — ► — oo. 
However, in the case where the cubic term is present, the solution u having the mirror 
symmetry about the y axis can exist because of the nonlinear effect. It will be explicitly 
shown in Sec. IV. 

III. NUMERICAL CALCULATION 

In this section, we explain how to calculate the solitonlike solution u of the cubic GS 
equation ( fl4l) . The numerical method is explained in Sec. IIIA. In order to establish the 
existence of the solitonlike solutions, we have to check the numerical errors carefully. This 
is discussed in Sec. IIIB. 

A. Numerical method 

In the numerical calculation, it is very convenient to choose the coordinates (X, Y) that 
are normalized by a quarter of the period in the y direction. For this reason, we introduce 
a parameter 

q := (tt/2)<z (16) 
and perform the coordinate transformation 

X:=q- X x, Y:=q- X y. (17) 

In the coordinates (X,Y), Eq. (HM becomes 

U t xX + U tYY = ~q 2 u(l + U 2 ), (18) 

and the period in the Y direction is Yp = 4. By the symmetries of the solution that we 
required in Sec. II, it is sufficient to solve in the range < X < X max and < Y < 1. 
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Here, X = X max is the outer boundary of the region of numerical calculation, and we choose 
X max = 5. The error coming from this cutoff value will be estimated in the next subsection. 

The boundary conditions are u = at Y = 0, uy — at Y — 1, and u t x = at 
X = 0. At the outer boundary X = X max , we have to impose the condition (fToT) . which is 
u = Aexp ^— X^Jix 2 /A — q 2 ^j sin(7r/2)y in the (X, Y) coordinates. Because we do not know 
the value of A before generating the solution of u, we calculate u y x and eliminate A. This 
leads to the so-called Robin boundary condition 



u x = —u 



(19) 



This formula is used as the boundary condition at X = X max . 

In order to solve Eq. ( fl8l) numerically, we adopted the second-order finite difference 
scheme with uniform grids. Since Eq. (j!8p is an elliptic equation, we have to prepare an 
initial surface and make it converge to the solution by the method of iteration. However, 
the solution was found to be unstable against this process. This is in contrast to the case 



of Ref. 



12| . where one of us solved a similar equation with no problem. The reason for 



the difference between the two cases is as follows. In both cases, the equation has the form 



V 2 u = —F(u). The function F(u) is monotonically decreasing as u grows in the case of 12), 
while it is a monotonically increasing function in the present found from Eq. (|18[) . 

Let us consider what happens in the latter case. Suppose the initial surface Uq is slightly 
larger than the real solution u. The program makes the surface approach the solution of 
the equation V 2 u = —F(uq). Because F(uq) > F(u), the solution of this equation u\ is 
further larger than uq, i.e. u\ > uq > u. Therefore, by continuing these processes, the value 
of u becomes larger and larger and eventually diverges. On the other hand, if u$ is a little 
smaller than u, the value of u becomes smaller and smaller and collapses to u = 0, i.e. the 
trivial solution. For this reason, the nontrivial solution of Eq. (|18|) is numerically unstable, 
and a new idea is needed to obtain it. 

Although we could not develop a new code which can automatically generate such nu- 
merically unstable solutions, we found a procedure to realize the sufficient convergence of 
iterations. The point is that for some initial surface uq, the surface u approaches the real 
solution u to some extent and then leaves it after that in the process of iteration. In other 
words, the real solution u behaves like an intermediate attractor in this process. Such a 
behavior typically occurs when uq crosses the real solution u, i.e. the regions uq > u and 
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uq < u both exist. Using this property, we proceeded as follows. We prepare a good initial 
surface and start the computation (the first trial). While the program is running, we 
observe the convergence parameters 



J2\ AU (IJ)\ A 

Ai := -, A 2 := max 



A«(7,J) 



(20) 



U(I,J) 

where (J, J) are the label of the grids and Auri t j) denotes the difference from the finite 
difference equation. The values of Ai^ first decrease and then increase. Just before Ai 
starts to increase, we write down the surface u^p and stop the program. Then, we prepare 
the new initial surface by = (1 + e)Up and run the program again (i.e. the second 
trial). Choosing e properly, we can make A 1]2 further smaller although a little experience is 
required in order to find the effective value of e. We continued these processes of trials until 
the conditions Ai < 10~ 8 and A 2 < 10~ 7 are achieved. 



B. Error estimates 



Using the above technique, we solved the cubic GS equation ( fl8i) for q = 0.1-0.9 with 
0.1 intervals. In all cases, we adopted the grid numbers (250 x 50). Since the solution is 
numerically unstable, we have to check the numerical error carefully in order to prove that 
our solution is not a numerical artifact. There are three sources of the numerical errors: 
the finiteness of X max , the finiteness of the grid sizes, and the truncation of the convergence 
process. 

The error by the finite X max value is evaluated by 5\ = [u(X max , l)] 2 , since the boundary 
condition ( fl5l) is derived by ignoring the cubic term in Eq. (TT4T) . The value of 8\ is less than 
0.02% for 0.1 < q < 0.7. It becomes larger as q is increased and 0.3% for q = 0.9. This is 
because the value of u decays very slowly for q ~ 1 by the boundary condition (floT) . We also 
compared the results of X max = 5 and 10 in the case q = 0.9. The error estimated in this 
way is 0.1%. 

The error by the finite grid sizes is estimated as 0.03% for all values of q by comparing 
the results of (250 x 50) and (125 x 25) grid numbers. This is natural because we used 
the second-order accuracy scheme and thus the error is expected to have the order of the 
squared grid size ~ 0.04%. 

The error by the truncation of the convergence process is estimated as follows. Suppose up 
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FIG. 1: 3D plot of the generated solution u for q = 0.5 in the region < X < 5 and < Y < 1. 

is the obtained solution and consider the equation V 2 w = — q 2 up(l + u 2 F ). If the convergence 
is perfect, the solution of this equation is u = up. Since the convergence process is truncated 
by the criterion explained above, the actual solution of u is different from up, and this 
difference indicates the error amount. In this way, the error is estimated to be less than 
0.03% for all values of q. 

Therefore, all the numerical errors are small and our results are reliable. 

IV. NUMERICAL RESULTS 

Now we show the numerical results. Figure Q] shows the 3D plot of the numerical solution 
of u in the range < X < 5 and 0<F<lforg = 0.5. The solution u(X, Y) takes its 
maximum value u pea k at (X, Y) = (0,1). Figure [2] shows the behavior of u(X, 1) on the 
line Y = 1 and Fig. [3] shows the behavior of u(0, Y) on the line X = (i.e. F-axis) for 
q = 0.1-0.9. The peak value u pea k increases as q is decreased. This is because the right hand 
side of Eq. ( TTBl) is proportional to q 2 and thus larger value of u is necessary for smaller q in 
order that the effect of the nonlinear term becomes relevant. From the right plot of Fig. [2l 
we see that the value of u(X, 1) decays more slowly for larger q because of the boundary 
condition (FKI1) . 

Figure H] shows the dependence of the peak value u pea k on q. From this figure, it is 
understood that u pea k diverges in the limit q — > 0. By plotting the relation between qu pea k 
and q, we found that u pea k is approximated by u pea k — 1.72/ q for small q. The solution of u 
becomes u = in the limit q — > 1, because u depends only on Y in this limit by the boundary 



9 





FIG. 3: The behavior of u(0,Y) on the F-axis for q = 0.1-0.4 (left) and q = 0.5-0.9 (right). 

condition (jl~9l) while we are solving the sequence for which u = at X = oo. By plotting 
the values of u^ eak as a function of 1 — q, we found that the formula u pea k — 2.27^1 — q 
approximately holds in the neighborhood of q = 1. 

We summarize the general properties that do not depend on specific forms of / and g. 
From Eqs. (TT21) and ffTTl) . the coordinates (X, Y) and (x, y) are related as X = a^q^x and 

20 
17.5 

15 
12.5 

Upeak 10 

7.5 
5 
2.5 



FIG. 4: The dependence of the peak value u pea k on q. The numerical data is shown by squares 
(□). u pea k behaves as u pea k cx q~ l and ^1 — q in the neighborhood of q = and 1, respectively. 
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Y = aoq l y. In the (x, y) coordinates, the period of \1/ in the y direction is 



yp 



Aq 

a 



(21) 



By the dimensional analysis, /(^) and g(^) are found to be expressed as 



Here, by Eq. (ED]), f(u) and g{u) are related as 



a 



o 



(22) 



f + 



u 



u z + — + C, 



where C is a non-negative constant. Calculating the magnetic field 
current (@J using Eqs. (ED, ( fl3l) and ( JTU1) , we obtain 

1 



A) 



- (-u x e x + u x e y ) + f(u)e z 



T 

A) 



— {-uye x + u, x e y ) + u(l + u 2 )e z 

q 



(23) 

and the electric 

(24) 
(25) 



From these formulas, the meanings of the parameters q, a , and /3 are understood. Since 
the inside of the parenthesis of Eq. ( |24"|) depends only on u and q (for a fixed form of /), 
the direction of the magnetic field at a given position (X, Y) is determined once the value 
of q is specified. This means that the shape similarity of the field lines is preserved when «o 
and (3q are varied. Furthermore, nondimensional quantities such as the beta ratio 2p/B 2 are 
independent of «o and (3q. Therefore, q is the parameter that determines all nondimensional 
properties of the system. For a fixed q, the value of a® determines the characteristic scale 
of the system through Eq. ( 12T1) . After fixing q and ceo, the magnitude of B is determined by 
specifying (3q. Hence (3q is (say) the field strength parameter. 

From x and y components of Eqs. (124H and (1251) . the magnetic field lines and the electric 
currents are confined on the contour surfaces of u. Figure shows the contours of u on the 
(X, y)-plane for q = 0.5. The directions of the magnetic fields are also shown. The magnetic 
fields are clockwise in the region u > and counter-clockwise in the region u < 0. From 
Eq. ( 1251) . it is seen that J z < in the region u > and J z > in the region u < 0. This 
relation between the directions of (B x , B y ) and the sign of J z is consistent with Ampere's 
law. 
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FIG. 5: The contours of |u| = 0.0-2.5 (with 0.5 intervals) for q = 0.5 on the (X, K)-plane. The 
arrows indicate the directions of magnetic field lines. They are clockwise in the region u > and 
counter-clockwise in the region u < 0. 

The z component of the magnetic field is specified by the function f(u). Changing f(u) 
affects (J x , J y ) through Ampere's law. If /(0) = and f(u) is a monotonic function in 
each region of u > and u < 0, a simple relation exists between the sign of B z and the 
directions of (J x , J y ). Let us consider the region u > 0. If f(u) is a monotonically increasing 
function, we have / > and f )U > 0. This indicates that 5 Z > and the electric currents 
are counter-clockwise. On the other hand, if f(u) is a monotonically decreasing function, 
we have / < and / )U < 0, which means that B z < and the electric currents are clockwise. 
The same relation is obtained also for u < 0. 



V. SUMMARY AND DISCUSSION 



In this paper, we studied the solitonlike solutions of magnetostatic equilibria by numeri- 
cally solving the cubic GS equation. Although the solutions were unstable against the nu- 
merical iteration, we found the procedure to realize the sufficient convergence and obtained 
the highly accurate solutions. The generated solutions are solitonlike in the x direction, 
periodic in the y direction and symmetric in the z direction. Our resultproves the existence 
of the solitonlike solution that was questioned in recent years j^, 3, 4, 5, (|. 

The solitonlike solution obtained in this paper behaves as an even function on a y = const. 
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line and has one extreme at the center. It is interesting to examine the existence of another 
solution that behaves as an odd function on a y = const, line and has an extreme in each 
region of x > and x < 0. Such a (say) 2-solitonlike solution could be expected by the 
following discussion. Denoting the obtained solitonlike solution by ^^'(x, y), the function 

¥ 2) ( Xl y) = - d/2,y) + d/2,y) (26) 

also approximately satisfies the GS equation (flOj) for sufficiently large d, since ^W(x) de- 
cays exponentially for large \x\. Therefore, one might expect the existence of 2-solitonlike 
solutions also for finite values of d. However, it is possible to show that no 2-solitonlike 
solution exists under the boundary condition \I/(0, y) = for any / and g satisfying 
/(0)f(0) + g'(0) = 0. To show this, we multiply to the GS equation © as 

# A + V,yyV, X = -(ff + Sf)**, (27) 

and integrate this equation over the region < x < oo and < y < yp. Assuming the 
exponential decay of \I/ at x 3> 1, the integrals of the second term on the left hand side and 
the right hand side vanish, and we have 

"VP 

* 2 J0,y)dy = 0, (28) 



o 

and therefore ty }X (0,y) = 0. Then, the GS equation ([9]) indicates that all derivatives of $ 
with respect to x vanish on the symmetry axis (assuming \1/ to be analytic). Hence \I/ = 
is the only solution. Physically, this means that when two or more solitons coexist, they 
interact each other and cannot be in equilibrium. 

It is interesting to discuss the stability of the solitonlike solution obtained in this paper. 
The system is expected to be unstable, since the magnetic islands are periodically located 
in the y direction and interactions between them are present. The most important factor 
for such interactions is the directions of the electric currents of the islands. If the currents 
of the islands are parallel (i.e., J z has the same sign), the islands attract each other and 

nn 

coalesce into larger islands. Such instability is known as the coalescence instability [13J, H4| . 
On the other hand, if the currents of the neighboring islands are anti-parallel (i.e., J z has 
an alternating sign), their interaction is repulsive and they tend to repel each other in the 
x direction as a result of small disturbance. Since J z has an alternating sign in our system 
as seen from Eq. ( 1251) . the repulsive instability is expected. In fact, both instabilities were 
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confirmed by the recent numerical work on the dynamics of magnetic islands with parallel 
and anti-parallel currents 6]. 

Although the plane-symmetric solitonlike solution in this paper could be of use in the 
contexts of the astrophysics or the solar physics, it would be more interesting to apply 
our method to the axisymmetric case. In the observations of active galactic nuclei (AGN), 
the astrophysical jets are often found to have knotty structures that suggest the presence 
of magnetic multiple islands j^]. Several models of the knotty jets have been proposed, 
and one of the possible directions is to model the knotty jets as magnetostatic equilibria 
0, i], [9]. Although these studies do not give the mechanism for the formation of the 
knotty jets, such models are expected to explain the long lifetime of collimation and knotty 
structure simultaneously. Namely, the knotty jets can maintain their shapes because they 
are in equilibrium in the comoving frame, and the time scale of the instabilities gives the 
lifetime of the knotty structure. The authors of (J studied the growth of instabilities of 
plane-symmetric solitonlike configurations by performing numerical simulations. Assuming 
that the plane-symmetric solitonlike systems well approximate the axisymmetric ones, they 
compared the results with the observations of the knotty jet of the radio galaxy 3C 303 
15, Q, [ljj]. Their conclusion is that the numerical simulation gives good agreement with 
the actual observations. Here, it should be pointed out that the assumption in that paper is 
not obvious and has to be justified. For this reason, the extension to the axisymmetric cases 
is necessary in order to examine if the solitonlike solutions can really model the astrophysical 
knotty jets. An axisymmetric quasiperiodic magnetostatic solution of the linear GS equation 
was proposed as the astrophysica. jet model fl. The numerical method in this paper enab.es 
us to generalize the study of [9( to the case of the nonlinear axisymmetric GS equation and 
thus to obtain further large class of astrophysical jet models as magnetostatic equilibria. 
The present work is the first step toward this direction, and we are planning to generalize 
our result to the axisymmetric case. It would be also interesting to further explore the 
solitonlike solutions in the helically symmetric cases. 
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